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INTRODUCTION 

Spectral collocation methods have proven to be efficient discretization 
schemes for many aerodynamic (see e.g., refs. 1-5) and fluid mechanic (e.g., 
refs. 6-9) problems. The high-order accuracy and resolution shown by these 
methods allows one to obtain engineering-accuracy solutions on coarse meshes, or 
alternatively, to obtain solutions with very small error. One drawback to these 
techniques has been the requirement that a complicated physical domain must map 
onto a simple computational domain for discretization. This mapping must be 
smooth if the high-order accuracy and exponential convergence rates associated 
with spectral methods are to be preserved (ref. 2). Additionally even smooth 
stretching transformations can decrease the accuracy of a spectral method, if 
the stretching is severe (ref. 5). A further difficulty with spectral methods 
has been in their implementation on parallel processing computers, where effi- 
cient spectral algorithms have been lacking. 

The above restrictions are overcome in the present method by splitting the 
domain into regions, each of which preserve the advantages of spectral colloca- 
tion, and allow the ratio of the mesh spacings between regions to be several or- 
ders of magnitude higher than allowable in a single domain. Such stretchings 
would be required to resolve the thin viscous region in an external aerodynamic 
problem. Adjoining regions are interfaced by enforcing a global flux balance 
which preserves high-order continuity of the solution, regardless of the type 
(diffusion- or advection-dominated) of the equations being solved. This inter- 
face technique maintains spectral accuracy, even when mappings and/or domain si- 
zes are radically different across the interface, provided that the discretiza- 
tion in each individual subdomain adequately resolve the solution there. Addi- 
tionally, the present technique allows spectral collocation methods to effi- 
ciently utilize parallel processing (ref. 10) where the application of conven- 
tional single-domain spectral discretizations have not been found to be effi- 
cient. 

A number of other spectral multi -domain techniques have appeared in the lit- 
erature. Application of finite element methodology, using Galerkin spectral 
discretization in the variational formulation within the elements, is a popular 
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technique (refs. 11-12). One drawback with such techniques is that a split 
Galerkin-collocation discretization must be used for convection-diffusion 
problems. These spectral -element methods also seem to work best in practice 
when used in a manner similar to classical finite element techniques: low-order 

internal discretization using many elements with no internal stretchings to 
improve resolution. In other words, resolution is increased in such methods by 
packing more elements in appropriate areas, rather than by increasing the order 
of discretization or applying an appropriate stretching within an element. The 
use of low-order spectral discretization over many elements does not take 
advantage of the exponential order convergence properties of spectral methods 
even if the split Galerkin-collocation discretization can display such 
properties. Techniques which the authors are aware of for interfacing 
collocation-discretized domains seem to involve explicit enforcement of contin- 
uity of the solution and its first derivative across the interface (refs. 

13-14). It is not clear how well these techniques perform for strongly convec- 
tion-dominated problems; the second author's experience with such techniques 
(ref. 1) has shown them to be not entirely satisfactory. 

Results 

In this section, the present global flux-balance spectral multi-domain method 
will be shown as applied to a number of one- and two-dimensional test problems. 
(Detailed explanation concerning the implementation of the algorithm is given in 
ref. 10.) The one-dimensional examples will serve to show that this method can 
maintain the exponential-order error convergence which is characteristic of 
collocation methods, even when adjoining domains have radically different dis- 
cretizations in terms of domain size, number of points, or stretchings. 

The first example is the equation 

U xx = cos Xe[ - 2 ’ 2] » u(_2) = u(2) = 0 (1) 

for which the exact solution is: 

U ( x ) = cos^ 5 (2) 

IT 

Eq. (1) is discretized in two unstretched domains: x^e [-2,0]; x^e 

[0,2], with N1 and N2 points, respectively. In table I is listed the maximum 
relative error as a function of N1 and N2. Note that the error decay is clearly 
faster than algebraic down to machine zero when N1 = N2. Also shown in table I 
is the behavior of the solution when one discretization is held fixed (Nl) and 
the other refined (N2). The overall error remains constant at a level essen- 
tially one-half of that seen when both domains were discretized at the coarser 
level. This is expected behavior, since the error in the interface condition 
has components from both domains; if one discretization is very much coarser 
than the other, its error will dominate the overall error of the solution. When 
both have equal error, then they must contribute equally to the overall error. 

The next example will illustrate the capability of the method for resolving 
very high gradients in a solution while imposing an interface condition which 
preserves spectral accuracy. Consider the viscous Burger's equation: 

u t + \ (U 2 ) x = v U xx ; x e [-1,1] (3) 

U ( -1 , t ) = U ( 1 , t ) = 0 U(x,0) = -sinUx) 

This problem has been studied extensively by a number of authors, using tech- 
niques ranging from standard finite difference, to single-domain spectral collo- 
cation and spectral element (ref. 15). The solution to this problem develops a 
very steep gradient region in the center of the domain; the slope at x = 0 
reaches a maximum, then decreases as the initial energy is dissipated away. For 
the parameters studied in ref. 15 (v = 0.01 /tt), this maximum is reached at t = 
0.5; a very accurate analytical solution gives a value of 152.00516 for this 
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maximum slope. The evolution of this solution calculated from the present 
method is shown in fig. 1 at time increments of 0.1. 

In the present study of this problem, three domains were used, the middle 
domain spanning a very small region (+0.05) around the "shock." Additionally, a 
mapping was applied in the middle domain to improve resolution. The maximum 
stretching allowable in this mapping is subject to the same restriction as 
stretchings in single-domain discretizations; e.g., maximum metric ratios on the 
order of 1(P. 

From the comparison study contained in ref. 15 the two methods giving the 
best accuracy for a given number of grid points were single domain spectral 
collocation and spectral element. The collocation scheme used a mesh stretching 
with a maximum-to-minimum metric ratio of about 100. Beyond this stretching a 
degradation in accuracy was found to occur. The spectral element discretization 
utilized four elements with 16 nodes in each. The behavior of the error in max- 
imum slope from these methods and the present scheme are shown in table II. As 
can be seen, the present method with just 35 total points (12 points in the 
outer domain, 13 points in the middle domain, 12 points in the left outer 
domain, hereafter denoted 12/13/12) yield results of equivalent accuracy to the 
spectral element and single-domain spectral collocation methods of ref. 16 both 
using 64 total points. Further mesh refinements using the present method show 
exponential-order error convergence, as seen in table II by the order-of- 
magnitude decrease in relative error as the mesh is refined to 20/21/20, and 
again with mesh refinement to 32/33/32. For the same total number of points, 
the present method is an order of magnitude more accurate than the single-domain 
collocation or spectral element solutions of ref. 15. 

In order to demonstrate the capability of the present method to handle 
radically-di fferent mappings between adjacent domains, a solution to the above 
viscous Burger's equation for v = 10" 4 , is shown in fig. 2. The maximum slope 
for this solution is greater than 5000. The discretization used was 12/31/12; 
the stretching in the middle domain was so severe that the ratio of largest mesh 
spacing in the outer domains to the smallest in the middle domain is greater 
than 10^. A factor of 5000 magnification of the high-gradient region of this 
solution is shown in fig. 3. The emphasis in this plot is the oscillation-free 
resolution of this region. (Linear interpolation between points is used for 
plotting, making the plot appear somewhat jagged.) 

To demonstrate the maintenance of conservation by the present interface 
technique, an initial condition was applied to the viscous Burger's equation to 
generate a moving "shock" which passed through the interface, as shown in fig. 

4. No oscillations, reflections, or abrupt changes in wave speed are seen as 
the "shock" passes through the interface. A very skewed discretization of 
12/17/27 was used for this case. Note that a multi -domain method formulated 
only for hyperbolic or for elliptic equations would be unable to perform well on 
this problem, since the dissipation-dominated region passes through the 
interface. 

Two-dimensional examples show similar performance of the present method. 

Shown in fig. 5 are constant-value contours for the solution of the Poisson 
equation: 

U xx + Uyy = cos^- cos^; x e [-2,2], y e [-1,1] (4) 

U ( x , — 1 ) = U(x,l ) = U(-2,y ) = U(2,y) = 0 

using two domains interfaced at x = 0 (the dashed line in fig. 5). Note that 
the contours pass smoothly through the interface. 

Table III contains the maximum relative error for various multi -domain 
discretizations of eq. (4). Exponential -order error convergence is again 
apparent from the rapid decrease with mesh refinements. The first group of re- 
sults in table III is for equal discretization in both domains, whereas the 
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second group shows the effect of having different discretizations in the 
direction along the interface. Only a small amount of accuracy is lost through 
the interpolation between two radically-different discretizations across the 
interface. 

An example of discretization on four domains with a corner point is shown in 
fig. 6, where isolines of the spectral solution to 

v 2 U = cos^ • cos^, x £ [-2,2], y e [-2,2] (5) 

U { x , -2 ) = U ( x , 2) = U(-2,y) = U(2,y) = 0 

are displayed. An influence-coefficient matrix algorithm described in ref. 10 
was used to compute this solution. This algorithm was developed for parallel- 
processor applications. 

Spectral discretization of problems with discontinuous coefficients or source 
terms (or equivalently, discontinuous transformation metrics), or discontinuities 
in boundary condition typically yields solutions with large oscillations and 
low-order error convergence. The present multi -domain technique may be used to 
isolate such discontinuities and recover exponential-order convergence. Figs. 7 
and 8 illustrate such an application; shown are solutions to Laplace's equation 
in which a jump in boundary conditions is enforced on one side of the domain. 
When the discontinuity lies at a point interior to one of the discretizations 
(fig. 7), oscillations are seen clearly in the solution isolines. When the dis- 
continuity occurs where the interface meets the boundary, however, the contour 
lines are smooth (fig. 8). 

Another example of the application of the multi -domain technique to isolate a 
discontinuity is in solution of the following equation: 

V • (kvU) =0; x e [-2,2], y e [-1,1] (6) 

U ( x , -1 ) = U(x,l) = U(-2,y) = 0 U(2,y) = cos ^ 

where k = kj, -2 _< x £ 0, and k = kg = 10 kj, 0 < x _< 2, with the inter- 
face at the line of coefficient discontinuity as shown in fig. 9. The computed 
solution is everywhere smooth, and the gradient jump at x = 0 is automatically 
enforced. 

To demonstrate the generality of the technique, eq. (4) was solved on the 
skewed two-domain mesh shown in fig. 10. This mesh, containing 17 x 16 and 18 x 
17 points in the left and right domains, respectively, was generated by first 
choosing the interface line, in this case a cubic polynomial. Chebyshev distri- 
butions with respect to arc length were used to establish the mesh points on the 
interface, as well as along the domain boundaries at x = +2. One curvilinear 
coordinate family was generated by connecting these corresponding points with 
straight lines. Mesh points along these coordinate lines were then established 
with Chebyshev distributions with respect to arc length, resulting in a sheared 
non-orthogonal mesh. Eq. (4) was written in generalized contravariant flux form; 
the metrics were evaluated by spectral differentiation of the coordinate distri- 
butions. The flux component normal to the interface was taken to be continuous 
in the interface condition. As can be seen in the isolines of the solution 
shown in fig. 11, the solution is everywhere smooth and regular. 

Conclusions 

The present global flux balance spectral multi -domain method has demonstrated 
maintenance of exponential -order accuracy herein on a variety of advection- and 
diffusion-dominated test problems. Extremely large difference in discretization 
across an interface, through domain size, number of points and stretchings, have 
been shown to not disrupt this property of the present method. Additionally, 
this technique can be used to isolate certain types of coefficient, mapping or 
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boundary condition discontinuities. A further application of the present method 
lies in implementation of spectral methods on parallel -processing computers, 
where the global nature of collocation methods have limited their effective- 
ness. The solution algorithm described is well -suited for machines with only 
nearest-nei ghbor connections between processing units. 

Further areas of applications being examined for the present method include 
implementation (1) in a time-dependent incompressible Navier-Stokes code for 
transition simulation, which will be run on the large scale Navier-Stokes 
Computer parallel computer under construction at Princeton University {private 
communication with D. Nosenchuck, Princeton University), and (2) in an 
external -aerodynamic compressible Navier-Stokes scheme which will interface a 
fine discretization of a nearfield region and a coarse discretization of the 
farfield to reduce storage and increase convergence rate. 
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Table I.- Maximum relative error vs. domain discretization for equation (1). 
Numerials in parentheses are powers of ten. 



5 , 5 

1.78 

( - 3 ) 

7 , 7 

9.81 

( - 5 ) 

9 , 9 

2.79 

( - 8 ) 

11,11 

4.85 

(- 11 ) 

13,13 

2.10 

(- 13 ) 

5 , 7 

8.82 

( - 4 ) 

5 , 9 

8.88 

( - 4 ) 

5,11 

8.88 

( - 4 ) 

7 , 9 

4.89 

( ~ 5 ) 

7,11 

4.90 

( - 5 ) 

7,13 

4.90 

( - 5 ) 

9,11 

1.39 

( - 8 ) 

9,13 

1.39 

( - 8 ) 
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Table II.- Maximum slope and percent relative error in maximum slope for viscous 
Burgers’ equation (eq. (3) ); comparison of present method with results from 
reference 15 • 


Method 

Discretization 

Max. slope 

% Relative error 

Present 

3 domains : 
12/13/12 

152.03544 

1.99 

(-2) 

Present 

20/21/20 

152.00011 

3.23 

(-3) 

Present 

32/33/32 

152.00513 

2.14 

(-4) 

Spectral 

element 

4 elements: 
16/16/16/16 

152.04 

2.29 

(-2) 

(ref. 15) 





Spectral 
colloc. 
(ref. 15 ) 

1 domain : 
64 

152.025 

1.31 

(-2) 

exact 
(ref. 15 ) 


152.00516 
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Table III.- Maximum error vs. domain discretization for equation (4) 


(N *1.V' 



( 6, 6), ( 6, 6) 

0.16809 

( -3) 

( 8, 8), ( 8, 8) 

.13260 

( -5) 

(10,10), (10,10) 

.71081 

( -8) 

(12,12), (12,12) 

.27290 

(-10) 

(16,16), (16,16) 

.71054 

(-13) 

( 6, 6), ( 6,12) 

.44939 

( -3) 

( 8, 8), ( 8,16) 

.54130 

( -5) 

(10,10), (10,20) 

.33158 

( -7) 
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Figure 2. - Computed solution to equation (3) at time increments of 0.1: 
v = 10-4, discretization 12/31/12, interfaces at + 0.02. 



Figure 3. - Expansion of high gradient region of figure 2. 
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Figure 5. - Computed solution to equation (4); interface at dotted line. 



Solution of Poisson eqn on Four Domains 


u xx + u yy = cos(7rx/4)'COs(7iy/4) 
D: xe[-2,2], ye[-2,2], u=0 on dD 



Figure 6. - Computed solution to equation (5) on four domains; interfaces at 
dotted lines. 
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Figure 10. - Skewed two-domain mesh. 
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Figure 11. - Solution of equation (6) on mesh shown in fig. 10. 
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